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Recent pulsar timing observations by the Parkers Pulsar Timing Array and European Pulsar 
Timing Array teams obtained the constraint on the relic gravitational waves at the frequency /« = 
1/yr, which provides the opportunity to constrain H„, the Hubble parameter when these waves 
crossed the horizon during inflation. In this paper, we investigate this constraint by considering the 
general scenario for the early Universe: we assume that the effective (average) equation-of-state w 
before the big bang nucleosynthesis stage is a free parameter. In the standard hot big-bang scenario 
with w — 1/3, we find that the current PPTA result follows a bound H t < 1.15 x 10 -1 mpi, and the 
EPTA result follows H* < 6.92 x 10 mpi. We also find that these bounds become much tighter in 
the nonstandard scenarios with w > 1/3. When w = 1, the bounds become H„ < 5.89 x 10~ 3 mpi for 
, the current PPTA and H, < 3.39 x 10 _3 mpi for the current EPTA. In contrast, in the nonstandard 

scenario with w — 0, the bound becomes H, < 7.76mpi for the current PPTA. 
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A stochastic background of relic gravitational waves, generated during the early inflationary stage, is a necessity 
dictated by general relativity and quantum mechanics [l|-[3| ■ The relic gravitational waves have a wide range spreading 
spectra, and their amplitudes depend only on the Hubble parameter in the inflationary stage, when the waves crossed 
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I. INTRODUCTION 



the horizon, and the expansion history of Universe after the waves reentered the horizon. So their detection provides 
a direct way to study the physics in the early Universe in both stages, during and after the inflation. 

Recently, there have been several experimental efforts to constrain the amplitude of relic gravitational waves in 
the different frequencies. The current observations of cosmic microwave background (CMB) radiation by the WMAP 
satellite place an interesting bound on the so-called tensor-to-scalar ratio r < 0.20 0, @, which is equivalent to the 
\ constraint on the energy density Q gw (f) of relic gravitational waves at the lowest frequency range / ~ 10~ 17 Hz. Among 
■ various direct observations, LIGO S5 has also experimentally obtained so far the most stringent bound fi gw (/) < 
\ 6.9 x 10 -6 around / ~ 100Hz Q. In addition, there are two bounds on the integration J Q gv Af)d\n f ^ 1.5 x 10~ 5 , 
obtained by the Big Bang nucleosynthesis (BBN) observation Q and the CMB observation Q. These bounds have 
been used to constrain the Hubble parameter (or the potential density of inflaton) in the inflationary stage, when the 
corresponding waves crossed the horizon [1, d, Q . 

The timing studies on the millisecond pulsars provide a unique way to constrain the amplitude of gravitational 
waves in the frequency range / G (10~ 9 , 10~ 7 )Hz [ltj. Recently, the Parkers Pulsar Timing Array (PPTA) team 
and the European Pulsar Timing Array (EPTA) team have reported their observational results on the stochastic 
background of gravitational waves and given the upper limit of J7 gw (/) at the frequency / = 1/yr [ill, fl2l|. In this 
paper, we shall infer from these bounds the constraint on H, the Hubble parameter at the waves' horizon-crossing 
time during inflation. In the calculation, we have considered a general early cosmological model, i.e., we assume the 
effective (average) equation-of-state w before the BBN stage can be of any value, which includes a wide range of 
cosmological scenarios. The derived bound of H would limit various inflation models. 

II. RELIC GRAVITATIONAL WAVES IN THE STANDARD HOT BIG-BANG UNIVERSE 

Incorporating the perturbation to the spatially flat Friedmann-Robertson- Walker (FRW) spacetime, the metric is 

ds 2 = a 2 (?7) [dif - (Sij + h l3 dx l dx j )] , (1) 
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where a is the scale factor of the universe, and r\ is the conformal time, which relates to the cosmic time by adrj = dt. 
The perturbation of spacetime hij is a 3 x 3 symmetric matrix. The gravitational-wave field is the tensorial portion 
of hij, which is transverse-traceless dih^ = 0, S^hij = 0. 

Relic gravitational waves satisfy the linearized evolution equation [l[ : 

d^V^g^hij) = -ltor&nj. (2) 

The anisotropic portion 7r,j is the source term, which can be given by the relativistic free-streaming gas [IH and the 
scalar field in the preheating stage [HI]. However, it has been deeply discussed that the relativistic free-streaming 
gas can only affect the relic gravitational waves at the frequency range / £ (10 -16 , 10~ 10 )Hz, which could be 
detected by the future CMB observations [l5j] . The generation of stochastic background of gravitational waves in the 
preheating stage has also been deeply studied (see, for instance, [HI]), where the gravitational radiation was produced 
in interactions of classical waves created by resonant decay of a coherently oscillating field. However, it was found 
that the typical frequencies of this kind of gravitational waves are quite high, i.e. / > 10 4 Hz. Even if the model with 
low energy H ~ lOOGeV is considered, the gravitational waves are important only at the frequency range / ~ 1Hz 
[3], which could be detected by the future laser interferometer detectors. So, both effects cannot obviously influence 
the relic gravitational waves at the frequency / £ (10~ 9 , 10~ 7 )Hz. For these reasons, in this paper we shall ignore 
the contribution of the external sources. So the evolution of gravitational waves is only dependent on the scale factor 
and its time derivative. It is convenient to Fourier transform the equation as follows: 

M**) = / £ [hkM^e^ + ac] , (3) 

where c.c. stands for the complex conjugate term. The polarization tensors are symmetry, transverse-traceless 
k l 4f(k) = 0, 6 ij e-f(k) = 0, and satisfy the conditions e^ lj (k)e^\k) = 26 ss < and e\f(-k) = e\f(k). Since the 
relic gravitational waves we will consider are isotropy, and each polarization state is the same, we have denoted 
h^\r]) by where k = \k\ is the wavenumber of the gravitational waves, which relates to the frequency by 

k = 27r/. (The present scale factor is set o,q = 1). So Eq. ([2]) can be rewritten as 

h k " + 2-h k ' + k 2 h k = 0, (4) 
a 

where the prime indicates a conformal time derivative dj dr\. For a given wavenumber k and a given time 77, we can 
define the transfer function tf as 

t f (r],k) = h k (r])/h k (r]i), (5) 

where 77, is the initial conformal time. This transfer function can be obtained by solving the evolution equation (|4|). 
The strength of the gravitational waves is characterized by the gravitational-wave energy spectrum, 

n gw = p gw /p Q , (6) 

where p gw = 32 * G {hijh 1 ^}, the critical density is p a = |^y, and H Q = lOO/i • km • ,s _1 • Mpc^ 1 is the current Hubble 
constant. Using Eqs. © and ([S]), the energy density of gravitational waves can be written as fl(| 

f dk Pt{k)i){ m ,k) 
Psw ~J k 32nG ' ij 

where Pt{k) = ^r\hk(i]i)\ 2 is the so-called primordial power spectrum of relic gravitational waves. Thus, we derive 
that the current energy density of relic gravitational waves 

/P t (k) 
dink figw(fc), and Ogw(A;) = Y^2*/(/?o,fe) (8) 

where the dot indicates a cosmic time derivative d/dt. 

Now, let us discuss the terms Pt(k) and tf(ji , k) separately. The primordial power spectrum of relic gravitational 
waves is usually assumed to be power-law as follows: 



k 



P t (k) = A t (K)(-) . (9) 
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This is a generic prediction of a wide range of scenarios of the early Universe, including the inflation models. Here, 
we should mention that there might be deviations from power-law if we consider the relic gravitational waves in a 
fairly large wave number span. In this paper, as a conservative consideration, we assume this form is held only when 
k is very close to the pivot wavenumber fc* . In the above expression, n t is the spectral index when k — > fc» . (n t = 
corresponds to the scale-invariant power spectrum.) A t (k t ) is directly related the value of the Hubble parameter H 
at time when wavelengths corresponding to the wavenumber fc* crossed the horizon [l], 0, [TtJ , 



aJ' (fc») 



'7T TOpi 



(10) 

k,=a,H, 



where mpi = 1/VG is the Planck mass. 

Now, let us turn to the transfer function i/, defined in ([5]), which describes the evolution of gravitational waves 
in the expanding Universe. From Eq. Q, we find that this transfer function can be directly derived, so long as the 
scale factor as a function of time is given. Actually, the analytical or numerical forms of tf have been discussed by a 
number of authors (see, for instance, fl8U2l| y 

In this paper, we shall use the following analytical approximation for this transfer function. It has been known 
that, during the expansion of the Universe, the mode function hk{ifj of the gravitational waves behaves differently in 
two regions [H||. When waves are far outside the horizon, i.e. k <C aH, the amplitude of hk keeps constant, and when 
inside the horizon, i.e. k 3> aH, the amplitude is damping with the expansion of Universe, i.e., hk ex 1/0(77). In the 
standard hot big-bang cosmological model, we assume that the inflationary stage is followed by a radiation dominant 
stage, and then the matter dominant stage and the A dominant stage. In this scenario, by numerically integrating 
Eq. ((4]), one finds that the damping function if can be approximately described by the following form |22h25| 



if(V0, k) = Z^£l ^ A + L36 ( * ) + 2 .50( * )*, (11) 
krjo il A U k eq k eq 



where k eq — 0.073r2 m /i 2 Mpc _1 is the wavenumber corresponding the Hubble radius at the time that matter and 
radiation have equal energy density, and 770 = 1.41 x 10 4 Mpc is the present conformal time. The factor VL m /Sl\ 
encodes the damping effect due to the recent accelerating expansion of the Universe 0, H3- I n t ms damping factor, 
we have ignored the small effects of neutrino free-streaming jl3| and various phase transition [2~l| . 
We can define a new function 



T(k) = t f (r]o,k)/yfl2Hl (12) 

So, the current density of relic gravitational waves becomes fL w (fc) = P t (k)T 2 (k) . In this paper, we shall focus on 
the wavenumber k k eq . In this range, we have [22I I24I |25| prjj 

15 fn m \ 2 1 



and a parameterized form for the current density of relic gravitational waves 

For the wavenumber k = fc*, the value of f2 gw (fc») depends only on the value of H*. So, in this standard scenario, 
an observational bound on the Q gw (k*) corresponds to a bound on the Hubble parameter H*, which will be shown 
clearly in Sec. [IVJ 



III. DAMPING FACTOR IN THE GENERAL MODEL OF THE EARLY UNIVERSE 

Although, in the standard hot big-bang universe, a radiation dominant stage is always assumed after the inflationary 
stage, there is no observational evidence to show this is held before the BBN stage. Actually, this assumption can be 
violated in a number of cases, for example, the existence of the reheating stage [HI, or the existence of the cosmic 



phase transition [21[. So, in general, before the BBN stage, one can assume that the average equation-of-state of the 



Universe is w, and the scale factor satisfies a simple power-law form 

a(xn 1+p . (15) 
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The constant fj relates to w by f3 = (—3w + l)/(3w + 1). Obviously, when w = 1/3, i.e. f3 = 0, it returns to the 
standard model. However, if the Universe is dominated by the kinetic energy of inflaton, one has w = 1 and (3 = —1/2. 
On the other hand, for a matter dominated era, one has w = and (3 = 1. 

Now, let us discuss the evolution of relic gravitational waves in this general cosmological model. In principle, it can 
be done by directly solving Eq. ((4]). In this paper, in order to avoid the complicated numerical calculation, we give 
an approximate method as below. 

We consider the wave h k with the wavenumber fc, which crossed the horizon at a = a k and the corresponding 
Hubble parameter is H k - So one has k = akHk/ag. One knows that, when the waves are in the horizon, h k oc l/a(n), 
damping with the expansion of the Universe, and when the waves are out the horizon, hk = constant, keeping its 
initial value. So one can define a ratio, which accounts for the damping of the gravitational waves, 

hk(r]o) a k a k a b 

777 — r = — = . (16) 

hkiVi) a o a b ao 

where at, is the scale factor at the temperature of Universe being lMeV, i.e. the BBN stage. 

In the standard model, where /3 = in (|T5|) is assumed (i.e. w = 1/3, the radiation dominant stage), we have 

where H b is the Hubble parameter in the BBN stage. Taking into account the relation k = a k Hk/aQ, we obtain that 

fefcfa) _ Ob / a b H b \ ^ 
hk(Vi) a o \ a k J 

However, in the general case with /3 ^ 0, we assume h k crossed the horizon at a — a k and the corresponding Hubble 
parameter is H k . (Note that, in general a k ^ a k and H k ^ H k , but k = a k H k /ao is still satisfied.) From the equation 
in (fT5")) . it follows that 

H k /a b \^ 



H b \a k 



So, in this general case, we have 



1+/3 



hkjrjo) _ fife _ at / a b H b \ 
hk(Vi) ao a Q \ a a k J 

Comparing Eqs. (| and (|T5)) . we can define the damping faction D(k) as follows 



(19) 



D(k) = (^4) / (20) 

V ilk [r]i ) j general \ rik ytji ) j standard 

(a Ho) {~k~) (fcl) (21) 
Thus, in this general scenario, the current density of relic gravitational waves becomes 

^gw(fc) = Pt{k)T 2 {k)D 2 (k), (22) 

which satisfies fi gw (fc) oc k nt ~ 2 ^ when k is close to k„. Using the formulae in ©, (jTUJ) , (fT3")l , (j2"Tj) . and substituting 
the cosmological parameters (h = 0.702, T C mb = 0.276K, Q A = 0.725, fl m = 0.275, and z eq = 3454 )H, we get the 
following simple result 

log 10 tt gw (M = 1-25 - + 21og 10 (-^), (23) 

6W + 1 TOpi 



where fc* = 2i:f stl and /» = 1/yr is used. In Sec. IIV1 we shall compare this with the observational results. 
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FIG. 1: The upper limit of f2 gw (fc») as a function of the parameter a. The black solid line (i.e. line 1) is for current PPTA 2a 
result [Tl|. the blue solid line (i.e. line 2) is for current EPTA 2a result [12fl . the blue dashed line (i.e. line 3) is for current 
EPTA la result Q2], and the red dotted line (i.e. line 4) is for future PPTA 2a result [Tp]. 



IV. CONSTRAINT BY THE PULSAR TIMING OBSERVATIONS 



Pulsar timing observations provide a unique opportunity to study the gravitational waves at the frequency range 
/ £ (10~ 9 , 10~ 7 )Hz. In 2006, Jenet et al. have analyzed the PPTA data and archival Arecibo data for several 
millisecond pulsars. By focusing on the gravitational waves with the wavenumer fc* (where fc* = 2irf* and /* = 1/yr), 
and assuming the density of gravitational waves satisfies il gw (k) = k 2+2a at around k ~ k*, the authors obtained the 
2a upper limit on f2 gw (fe*) as a function of a [ll[, which has been shown in Fig. [1] (black solid line). This figure shows 
that f2 gw (fc*) < 4.05 x 10~ 8 when a = — 1. However, this upper bound increases to be 1.98 x 10 -6 when a = 0. 

Recently, this upper limit has been updated. In [lj], the authors have used the current data from the EPTA to 
determine an upper limit on the stochastic gravitational-wave background as a function of the spectral slope a. The 
la and 2a bounds are shown in Fig. [T](blue lines), which are slightly lower than those in PPTA case for any given a. 

It is interesting that in the authors have also investigated the possible upper limit (or a definitive detection) 
of stochastic background of gravitational waves by using the potential completed PPTA data-sets (20 pulsars with 
an rms timing residual of 100 ns over 5 years). We have also plotted this potential upper limit in Fig. Q](rcd dotted 
line) . 

Now, let us compare these observations with the analytical formulae of relic gravitational waves in Sec. IIIII Firstly, 
it is necessary to relate the parameter a with the theoretical models. In Sec. IIIII Eq. (|2"2"|) shows that il gw (k) oc fc Ilt ~ 2/3 , 
where fi = (— 3w+l)/(3w+l). Comparing this with the assumed form O gw (/c) ex k 2+2a , we get the interesting relation 

n t 2 

a = — . (24 

2 3u; + l y ' 

This relation shows that, in the standard hot big-bang scenario with w = 1/3, and the scale-invariant primordial 
power spectrum with n t — 0, we have a = — 1. In this case, let us use the bounds of gravitational waves to constrain 
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FIG. 2: The upper limit of the Hubble parameter H* as a function of the parameter a (left panel) and the average equation- 
of-state w (right panel), where we have assumed nt = 0. In both panels, the black solid lines (i.e. line 1) are for current PPTA 
2a result, the blue solid lines (i.e. line 2) are for current EPTA 2a result, the blue dashed lines (i.e. line 3) are for current 
EPTA la result, and the red dotted lines (i.e. line 4) are for future PPTA 2a result. 



the Hubble parameter fZ* in the inflationary stage. Taking into account the formula in Eq. (|2"3")l and using w = 1/3, 
we obtain the 2er upper limit of -H*, i.e. < 1.15 x 10 _1 mpi for the current PPTA case, iJ* < 6.92 x 10 _2 mpi for 
the current EPTA case, and the future PPTA is expected to give fz* < 7.94 x 10 _3 mpi. These results are listed in 
Table H 

Although, the inflation models always predict the nearly same Hubble parameter throughout the inflationary stage, 
it is necessary to constrain H , the Hubble parameter at quite different stages of inflation, which encodes the evolution 
information of inflaton. Here, let us compare the bound of H inferred from pulsar timing with those obtained in CMB 
observations and LIGO observations. The recent CMB observations by the WMAP satellite provide the constraint 
on the tcnsor-to-scalar ratio r < 0.20 HI, which is equivalent to the bound of H/mp\ < 6.92 x 10 -6 , where H is the 
Hubble parameter of inflation when the waves with frequency / = 1.94 x 10 _17 Hz crossed the horizon. The recent 
LIGO 5S reported so far the tightest constraint f2 gw (/) < 6.9 x 10 -6 on relic gravitational waves at the frequency 
/ ~ 100Hz J6[, which corresponds to H/nipi < 1.46. Comparing with these results, we find that the current and the 
potential future pulsar timing constraints on H are quite tighter than that of LIGO, but much looser than the CMB 
constraint. 

Now, let us relax the assumptions of the early Universe. We only assume n t = 0, which is approximately held in a 
wide range of inflation models. So, we can constrain the Hubble parameter in a wide range of w by the following 
inequality 

log 10 (— ) < \{U{a) - 6.74a - 1.25), (25) 
mpi 2 

where U(a) is the upper limit of f2 gw (fc ! „) based on the pulsar timing observations, which is a function of the parameter 
a (in this case, a relates to w by the relation a = — 2/(3w + 1)). The bounds of iJ* as functions of a (left panel) 
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FIG. 3: The upper limit of the Hubble parameter H* as a function of the average equation-of-state w, where we have considered 
three cases with nt = —0.5, 0, 0.5. The black lines are for current PPTA 2a result, and the red lines are for future PPTA 2a 
result. 



TABLE I: The 2a upper limit of the quantity H*/mp\ inferred from various pulsar timing observations. In this table, we have 
assumed rit — 0. 





Current PPTA 


Current EPTA 


Future PPTA 


w = 


7.76 




1.41 


w = 1/3 


1.15 x 10" 1 


6.92 x 10 -a 


7.94 x 10- 5 


w = 1 


5.89 x 10 _:5 


3.39 x 10 _a 


3.55 x 10" 4 


w — ► oo 


3.16 x 10~ 4 


1.20 x 10~ 4 


1.41 x 10~ b 



and w (right panel) are shown in Fig. [2] These bounds in three special cases with w = (i.e. a = —2), w = 1 (i.e. 
a = —1/2) and w — > oo (i.e. a = 0) are also listed in TablelH Clearly, we find that a larger w corresponds to a tighter 
bound of Especially, in the limit case with w —> oo, the current EPTA gives the constraint < 1.20 x 10 _4 mpi, 
and the future PPTA is expected to give a bound of < 1.41 x 10~ 5 mpi. 

In the end, let us discuss the most general case with free parameters n t and w. In this case, the inequality (|2"5)) 
becomes the constraint on the physical parameters n t and 7J* as follows 

TT 1 

1 ogin( — ) - l-69n t < -(U(a) - 6.74a - 1.25). (26) 
mpi 2 

Here, we should remember that a relates to the physical parameters by Eq. (|24j) . The spectra index n t influences 
the bound of mainly by slight changing the corresponding relation between a and w. In Fig. [31 we calculate the 
upper bound of in two special cases with n t = 0.5 and n t = —0.5, and compare them with those in the case of 
n t = 0. This figure shows that the parameter n t only slightly affects the bound of and a larger n t follows a looser 
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bound of iJ*. For example, the current PPTA observations follow < 1.82 x 10 1 mp\ at the case with n t = 0.5 and 
w = 0, which is only 1.6 times larger than the bound < 1.14 x 10 _1 mpi at the case with n t = and w = 0. 
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